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ABSTRACT 

Numerical simulations are becoming a more effective tool for conducting de- 
tailed investigations into the evolution of our universe. In this article, we show 
how the framework of numerical relativity can be used for studying cosmolog- 
ical models. The author is working to develop a large-scale simulation of the 
dynamical processes in the early universe. These take into account interactions 
of dark matter, scalar perturbations, gravitational waves, magnetic fields and a 
dynamic plasma. The code described in this report is a GRMHD code based 
on the Cactus framework and is structured to utilize one of several different dif- 
ferencing methods chosen at run-time. It is being developed and tested on the 
Texas Learning and Computation Center's Xanadu Cluster. 

Subject headings: Cosmology: Numerical Relativity: Rclativistic Plasma 



-3 - 



Introduction 



Our knowledge of how the universe evolved comes primarily from observations of large 
structures such as stars, galaxies, clusters and super-clusters of galaxies as well as from 
observations of the Cosmic Microwave Background (CMB) Radiation. Based on these 
observations, the Standard Model of Cosmology was developed during the mid to late 
twentieth century. Some elements of this model include the existence of primordial metric 
perturbations, m agnetic fields and an ea rly universe filled with a nearly homogenous and 



isotropic plasma (jMukhanov et al. 



19921 ). The Perturbed Friedmann- Robert son- Walker 



(FRW) metric, which describes the space-time curvature of the early universe, takes the 
following form, 



ds' 



[1 + 2(j))dt 2 + u i dtdx i + a{tf [((1 + 2 ip) + h ij )d^dx i ]. 



(1) 



Here a(t) is the scale factor and 0, if), Ui, and hij are the scalar, vector and tensor 
perturbation terms. Many cosmological models relate density fluctuations and variations 
in the CMB to perturbations in the FRW metric at the time of recombination. These 
perturbations start off small and grow as a power-law with time as the competing forces o f 



universal expansion and gravit ational attraction affect t heir growth (IMukhanov et al. 



Work by Kodama and Sasaki (IKodama & Sasaki 



19841) . Sachs and Wolfe flSachs fc Wolfe 



1992). 



1992J) all showed 



19671 ) and Mukhanov, Feldman and Brandenberger (jMukhanov et al. 
analytically how metric perturbations could cause density perturbations in a hydrodynamic 
fluid. 



2008aJ) 



Recently, beyond the Standard Model cosmological theories (jKahniashvili et al 
have suggested that primordial fluids and fields are also potential sources of observable 
gravitational waves. This realization opens up exciting new possibilities, making primordial 
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gravitational radiation an important source of information about the early universe. Taken 
together, this leads to the idea that there was a dynamical interaction between matter, 
electromagnetic and gravitational fields in the early universe that affected the evolution of 
the universe. Signatures of these interactions may still be observable today. The objective 
of this paper is to show how the tools of numerical relativity can be used to study such an 
interaction and to introduce a computer code written for this purpose. 

This project uses the framework of numerical relativity to develop a computational 
laboratory to study the evolution of the early universe. The initial focus of this work is 
on deriving the spectrum of gravitational waves produced by relativistic turbulence in the 
early universe. Future work may involve a more advanced study of how this gravitational 
wave spectrum is effected by the presence of dark matter or a pre-existing primordial 
gravitational wave field. Numerical Relativity has been used for years to study the collisions 
of compact objects such as black holes and neutron stars and to predict the spectrum of 
the gravitational waves produced by their interactions. We propose to use this tool to 
provide detailed studies of cosmological events that may also one-day be observable using 
either gravitational or conventional astronomy. In cosmology, numerical simulations are 
capable of providing more detail than the analytic calculations that have been performed 
to date. For example, a modern General Relativistic Magnetohydrodynamic (GRMHD) 
code is capable of performing such simulations by evolving both the spacetime and plasma 
field dynamically and can therefore represent chaotic processes such as turbulence more 
accurately than analytic calculations. This paper is a summary of the techniques used to 
develop such a simulation. 
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2. Development of Initial Conditions 

Every numerical simulation consists of three parts: Initial Data, Numerical Evolution 
and Data Analysis. In this section, we will focus on the initial data needed for cosmological 
studies. 



2.1. Background and Hypotheses 



Several different mechanisms for producing primordial gravitational waves have 
been identified and studied by researchers. These include quantum fluctuations during 
inflation, bubble wall motion and collisions during phase transitions, cosmological 
magnetic fields, osc illating classical fields during reheating, cosmological defects an d 



plasma turbulence ( 



Caprini fc Durrei 



2006 



Caprini et al. 



2009 



Dolgoy et al. 



2005bi 



2002 



Gogoberidze et al 



200 



Nicolis 



200, 



Kahniashvili 



2005a 



Kahniashvili et al. 



2006 



2008aJJb, 



20041 ). We assume that the early universe was in a metastable state during a 
first order phase transition. That false vacuum was separated from the true vacuum by a 
potential barrier or a scalar field. Quantum tunneling occurred across the barrier in finite 
regions of space resulting in true vacuum bubbles inside the false vacuum phase. As the 
universe expanded and cooled, the energy difference between the false vacuum and true 
vacuum got larger, making the phase transition more probable. Eventually, the probability 
of nucleating one critical bubble per Hubble time became high enough to cause the phase 
trans ition to begin. This define d the transition temperature, which is believed to be about 1 



TeV (jKahniashvili et al. 



2008af ). The nucleated bubbles expanded and collided, eventually 
filling the whole universe. The collision of two or more of these bubbles broke spherical 
symmetry and released some of their energy as gravitational waves. Since the expansion 
of the bubbles was accompanied by macroscopic motions in the cosmic matter field, the 
collision of these bubbles also resulted in the anisotropic stirring of the field. This caused 
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turbulent motions which provided a primary source of gravitational waves for this research. 



2.1.1. Primordial Magnetic Fields 



Magnetic fields are believed to have played a large part in the dynamics of the universe's 
evolution. Unfortunately, the only thing that we currently know about magnetic fields in the 
early universe is that little is known about their existence or absence. There are no direct 
observations of primordial magnetic fields. Our only knowledge of primordial magnetic 
fields comes from theories that fail to agree on their amplitude. There are currently 



several dozen theories abou t the origin of cosmic magnetic fields (IBattaner fc Lesch 



2000 



Grasso &: Rubinstein 



20011 ) . The main reason that we believe that primordial magnetic 



fields existed is because they may have been needed to seed the large magnetic fields 



observed t oday. Most theories of cosmic magnetic field generation fal 



categories ( IBattaner fc Lesch 



2000 



Dolgov 



2001 



Grasso Rubinstein 



into one of three 



20011 ): 1) magnetic 



fields generated by phase transitions; 2) electromagnetic perturbations expanded by 
inflation; and 3) turbulent magnetofluid resulting in charge and current asymmetries. 

Most models calculate the magnitude of primordial magnetic fields by starting with 
the observed strength of galactic or intergalactic magnetic fields and calculating how this 



field should have been amplifie d or diffused by external effects such as the ga 



and expansion of the universe (IBattaner Sz Lesch 



2000 



Grasso fc Rubinstein 



actic dynamo 



20011 ). A 



major problem is that there doesn't appear to be a universal agreement of how efficiently 
a galactic dynamo could have strengthened seed magnetic fields. Estimates of the strength 
of these seed fields can vary by tens of orders of magnitude. Seed magnetic fields produced 



during Inflation are predicted to have a current strength somewhere between 



10 9 G on a scale of a few Mpc ( IBattaner &: Lesch 



2000 



Grasso fc Rubinstein 



0- 11 G and 



2001 



Hoyle 



19511 ). Magnetic seed fields generated by phase transitions are believed to be less than 10 



-23 
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G at galactic scales (IBattaner &: Leschll2000l ; iGrasso &: Rubinstein! 120011 ). Some turbulence 



theories imply that magnetic fields were not 



therefore requiring no magnetic seed fields (IBattaner fc Lesch 



generated until a: 



ter the first stars were formed 



2000) 



Given how little is understood about primordial magnetic fields and the general lack 
of agreement among theoretical predictions, it seems clear that the existence of primordial 
magnetic fields can neither be confirmed or ruled out. It seems that the best we can 
do is set an upper limit on the strength of primordial magnetic fields and utilize this 
limit as a starting point in developing models of cosmic turbulence. Observations of the 

intensity of the magnetic seed fields to a 



current value of 10 9 G ( 


Battaner & Lesch 


2000 


Muravama & Pierce 


2002 


)• 



Grasso fc Rubinstein 



2001 



Hoyle 



1951 



It is well known that gravitational waves can i nteract with a mag netofluid in the 



2005b! ) showed that 



presence of a magnetic field. Work by Duez et al (IDuez et al. 
gravitational waves can induce oscillatory modes in a plasma field if magnetic fields 
are present. Work by Kahniashvili and others have shown how a turbulent plasma can 
yield gravitational waves. The result may be a highly nonlinear interaction as energy is 
transferred from the fluid to the gravitational waves and back. 



2.1.2. Theory of the Electron) eak Phase Transition 



Our treatment of early universe phase transitions is based on fundamental particle 
physics, rather than models. For the EWPT we introduce one supersymmetric field, the 
Minimal Supersymmetric Standard Model (MSSM), to obtain a first order phase transition. 
Finding the quantum of this field, the Stop, would be a major success at the LHC. 



Our electroweak theory (IHenley et al. 



20101 ) is based on the standard Weinberg-Salam 
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model with the addition of the Stop, the supersymmetric partner to the top quark. The 
Lagrangian is 

c mssm = £i + £2 + C 3 + i eptonic anc i quar ^ interactions (2) 



& = ~\w; u W^ - \b, v B^ (3) 



C? = \(id, - f A<^)$ s | 2 - V hs ($s, «&) , (5) 

where the W\ with i = (1,2), are the W + ,W~ fields, is an SU(3) gauge field, ($, $ s ) 
are the (Higgs, right-handed Stop fields), (r\X a ) are the (SU(2),SU(3)) generators, and the 
electromagnetic and Z fields are defined as 

A e ™= J— ig'W' + gB,) (6) 
V 9 +9 

Z,= 7 J== (gW*-g'B,). (7) 
Vsr + 9 2 

The equations of motion are obtained by minimizing the action using the magnetic 
fields created by bubble nucleation. The bubble collisions res ulting from the EWPT 



20091 ). 



nucleation were used in the prediction of gravitational waves ( iKahniashvili et al. 
This fundamental theory will be the basis of our research on primordial magnetic fields and 
gravitational waves associated with the EWPT. 



2.1.3. Theory of the Quantum Chromodynamic Phase Transition 

For the Quantum Chromodynamic phase transition (QCDPT), an essential aspect of 
our proposed program of creation of gravitational waves, we start from basic QCD, as the 
QCDPT has been shown in recent studies to be first order, with bubble collisions creating 
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magnetic fields. The interactio n Lagrangian us ed to derive the magnetic fields created by 



QCDPT bubble wall collisions ( iKisslinger 



20031 ) is 



C mt = -e^A e ™V, (8) 

where \I/ is the nucleon field operator. This leads to the electromagnetic interaction with 
the nucleons magnetic dipole moment given in terms of the electromagnetic field tensor, 
F^ v by 

V mt = 2^<V75*i^ . (9) 

For the gluonic instanton wall oriented in the x-y direction one obtains for 
B z = Bw = F 21 within the wall of thickness p 

A suppression factor of (pAqcd) -1 has been used in Eq. ffTUl) since the aligned dipoles tend 
to cancel. Therefore in our picture, at the end of the QCD chiral phase transition, there is 
a magnetic wall in the hadronic phase, which can be the seed for galacti c and extra-galactic 



magn etic fields. This was used to estimate gravitational wave creation ( jKahniashvili et al. 



20091 ). This is the basis for our proposed research on magnetic field evolution and 
gravitational wave creation in the hadronic phase, after 10 -6 seconds. 



2.1.4- Turbulence in the Early Universe 

Turbulence provides a particularly interesting GW source because it is not well 
understood analytically. This turbulence is a natural result of dynamics of the early 
universe resulting from bubble wall collisions and other chaotic events during the first order 
phase transitions. Analytic work done to date summarizes the dynamics of the first order 
phase transitions using two quantities, a and (3. a is traditionally defined as the ratio of 
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false vacuum energy and plasma thermal energy density This provides a measure of the 
transition strength. If a is much less than one, the transition is very weak. If a is larger 
than unity, the transition is very strongly first order. j3 is the rate of variation of the 
nucleation rate at the transition time. It fixes the time scale of the phase transition once 
the transition has begun. After a time interval f3~ 1 , the whole universe is converted to a 
true vacuum phase. Therefore the turbulent stirring should only last f3~ 1 . 

The amount of gravitational waves emitted by bubble collisions and turbulence 
generated in the plasma are also determined from two quantities, k and tv k is the fraction 
of vacuum energy transferred into fluid kinetic energy and v b is the velocity of bubble wall 
expansion. Bubble walls can propagate via two modes, detonation and deflagration. For 
detonation, the bubble walls are thin compared to the radius and they propagate faster 
than the speed of sound. This results in: 



1/V3+ {a 2 + 2a 



/3U/2 



v b (a) = IV '\, ^— (11) 

1 + a 



k(ch) = 



4 3a 
0.715c + -J Y 



(12) 



1 + 0.715a 

If the bubbles propagate by deflagration, the walls are thick and have a lower energy 
density. It is currently believed that for a relativistic plasma, the deflagration expansion 
mode is unstable so only the detonation modes will result. 

The number density of turbulent eddies within a Hubble radius should depend on Vb 
and f3. The characteristic velocity perturbation of the turbulent fluid for the largest eddies 
at the stirring scale is given by: 

/ 3na , . 

"° = VlTw (13) 

For kol ~ 1, which corresponds to a strongly first order phase transition, vq is about 0.65 at 
the time of the electroweak phase transition. 
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2.2. Computational Model 



In order to study the interaction of the plasma field and the background spacetime 
dynamically (or separately) our team has written and is testing/improving a General 
Relativistic Magneto-hydro dynamic (GRM HD) model implemented using the Cactus 



computer code framework ( jGarrison 



20081 ). We now describe the basic variables and 



equations that constitute this model. 



2.2.1. The Spacetime Evolution Model 



The spacetime metric can be written as: 



ds 2 = -N 2 dt 2 + 7i i (f, t){dx l + N l dt)(dx 3 + N ] dt). 



(14) 



Here N is the lapse, N l is the shift vector and 7^ is the spatial 3- metric (IBaumgarte fc Shapiro 



1999 



DuezetaL 



2005al ). For this work, 3-metric and its "time-derivative", the extrinsic 



curvature, will be evolved using a str ongly hyperbolic version of the BSSN formulation 



of numerical relativity ( IBrown et al 



2012|). 



2.2.2. The General Relativistic Magnetohydrodynamic Model 



The fluid and electromag netic fields of the GRM HD equations are developed from 



several well-known equations ( IGammie &: Toth 



20031 ) . They include the conservation 



of particle number, the continuity equation, the conservation of energy-momentum, the 
magnetic constraint equation and the magnetic induction equation. For a system consisting 
of a perfect fluid and an electromagnetic field, the ideal MHD stress-energy tensor is given 
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by 



b 2 

T" u = ( Po h + b 2 )u»u v + (P + — ) - U L b v (15) 
P 

h = l + e+ — (16) 
Po 

(17) 



If = 




B (u) = 


a 


B U = 


1 B l 

n v 



(«) 

(18) 

+ £(V) • ( 19 ) 

Here, P is the fluid pressure, po is density, B % is magnetic field, ti M is four- velocity, is the 
enthalpy, e is specific internal energy, and b 2 is the magnitude of the magnetic vector field 
squared. The addition of viscosity modifies the MHD stress-energy tensor by incorporating 
the viscous stress tensor 

b 2 

T» v = (p h + b 2 + Q)u»u v + (P + — ) g» v - W + S^. (20) 

Here Q is artificial bulk viscosity and Y^ v is the viscous stress tensor for artificial shear 
viscosity In order to keep our model as true to reality as p ossible we utilize art ificial 



2005J) are 



viscosity as a way of handling shocks. Our viscosity terms (lAnninos et al. 
described defined below and are only meant to be used when the divergence of the fluid 
flow is negative. 



Q = I n Al d k V k (k q Al d k V k - h C s ) (21) 

d k V k 
3 

I n = (p + pe + (P + Q + b 2 )) N. (23) 



E} = I n Al{k q Al d k V k - h C s )Sym(5 J V l - 5j) (22) 



In these equations V % is the fluid velocity, C s is the local speed of sound, Al is the minimum 
covariant zone length and k q and ki are constants multiplying the quadratic and linear 
contributions, respectively. The Sym(...) function in the shear viscosity equation is a 
symmetry operation. 
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Dark matter can be added to the system using a two-fluid approach where the 
stress energy tensor for dark matter is added directly to the stress energy tensor for the 
magnetofluid, therefore, completing the right-hand side of Einstein's equation. 

2.2.3. Initial conditions: Plasma Field 

In order to get a more accurate picture of the primordial gravitational wave spectrum, 
we directly simulate the turbulent primordial universe with the most realistic initial 
conditions possible. These should include not only a plasma field with a realistic EOS, 
they should also include elements such as magnetic fields, dark matter and possibly even 
gravitational waves produced by other sources. 

The study will begin at t > 1CT 6 seconds after the Big Bang near the beginning of the 
Hadron Epoch, when the primordial plasma field began to look like a relativistic plasma and 
the strong force could be safely ignored. At this point the Debye length is about 1CT 16 m so 
even a small computational domain should demonstrate the dynamics of the plasma. The 
plasma at this time was composed mainly of electrons, positrons, neutrinos and photons. 
Many of the initial conditions at this epoch are well known or can be fairly easily calculated 
using available literature. In addition, given that most cosmological models agree that over 
80% of the matter in the universe is composed of dark matter, our cosmological simulations 
can also include some non-magnetized "dark" fluid. We take this into account by adding a 
second pressureless non-magnetized fluid to our initial plasma field. 

The initial matter field will be taken to be homogenous and turbulent. We introduce 
turbulence into the system by randomly varying the initial velocities of the fluid elements 
up to the magnitude of vq. In addition, we will be working to better establish the initial 
conditions resulting from the first order EWPT. Based on the arguments in the previous 
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2000 



sections (IBattaner fe Lesch 



Maroto 



2000 



Dolgov 



2001 



Grasso fc Rubinstein 



2001 



Jedamzik et al. 



20011 ) . at this time initial magnetic field should be less than or equal to 10 
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G. 



2.2.4- Initial conditions: Spacetime 

For this computational study the initial spacetime is constructed in such a way as to 
mimic the conditions present during the Hadron epoch (t > 10~ 6 seconds). I choose to 
begin the simulation during the Hadron epoch because at that time the primordial plasma 
field appeared to look like a relativistic plasma field that could be modeled using a GRMHD 
code as opposed to a Quark Gluon plasma field. At this time the strong force could be 
safely ignored as electromagnetic effects would dominate the plasma's dynamic motions. 
Also, at this time any EWPT would be complete. 

Although initial calculations will use a fixed spacetime background, we may find it 
necessary to dynamically model the turbulent/spacetime interactions in order to fully 
understand the physics of the early universe. To do this, we need good initial conditions for 
the curvature of space during this epoch. 

The Robertson- Walker (R-W) spacetime metric for a flat universe can be written as 

ds 2 = -dt 2 + a 2 [t) (24) 

where t is the timelike coordinate, g^ is the maximally symmetric three-dimensional space 
metric, and &(t) the scale factor. For calculating the initial spacetime, conformal time, 
r = — , is often used instead of cosmic time, t = t a 2 , to simplify the equations. 

Therefore, 

ds 2 = a 2 (r) g^dxPdx" (25) 
Here, the scale factor and Hubble parameter can be calculated based on the temperature 
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and mass-energy density of the early universe relative to today In order to add gravitational 
waves, the metric g^ u is broken into two components: g^, the background metric plus a 
perturbation, h^ v . The metric can be written as: 

9iiv = g^iu + hpu. (26) 

The background metric may take any form. It is also assumed that any perturbations are 
linear. For this study, a R-W metric is presumed. Initial perturbations may or may not be 
used for numerical experiments within this study. This metric may involve scalar, ((f), ijj) 
vector (ui) and tensor (hij) perturbations. The tensor perturbations are symmetric and 
transverse-traceless so dih^ = and 8^ hij = 0. The initial amplitude and spectrum of any 
of these fluctuations depends on the theory used to explain the generation of perturbations 
from inflation. Note that by comparing equations (TH|) and ([I]), it can be shown that 
scalar (0) and vector (ui) perturbations can be related to the chosen lapse and shift in the 
same way that scalar (-0) and tensor (hij) perturbations can relate to the three-metric and 
extrinsic curvature. Because the focus here is on tensor perturbations, a Geodesic Slicing 
is used so the lapse, N, is set to unity and the shift vector, N l , is set to zero. Later, if 
scalar perturbations are included, they can be added by modifying the lapse, as well as the 
three-metric and extrinsic curvature. 

The process of generating tensor perturbations or gravitational waves from quantum 



fluctuations during inflation is similar to the pr ocess of generating sea 



pertu rbations. For example, work by Grishchuk (1 Grishchuk 



1998 



2001 



ar or vector 



Grishchuk et al. 



20011 ) gives a basis for calculating the spectrum and amplitude of these waves for different 



slow-roll parameters. 

The possibility of all polarizations are included using hfj, hfj, hfj and hfj. Here, hfj and 
h*j are the plus and cross polarizations respectively, and hfj and hfj are the left and right 
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rotating polarizations defined by: 



h^-Lihfj-ih*), h ^^=(h± + ih^. (27) 



According to work by Alexander ([Alexander &: Martini 120051 ) . rotating polarizations should 



dominate in the early universe and satisfy the equations: 

Uh% = -2ieh'^ Uhf 3 = +2i9h!*. (28) 
Here, the dot denotes a time-like derivative with respect to r and the prime denotes a spatial 



derivative along the gravitational wave's direction of propagation. I 



set to zero, the unpolarized gravitational wave signature is recovered (j Garrison fc de la Torre 



Alexander's 6 value is 



20071 ). For this work, rotational polarizations may have the added benefit of introducing 
extra vorticity into the homogeneous plasma. This may result in an increased magnetic 
field due to the dynamo effect. During inflation, a non-zero 6 term results in a decrease in 
the amplitude of h L and an amplification of h R and, therefore, cosmological birefringence. 

Initial gravitational wave spectra from sources such as bubble collisions or phase 
transitions can be added linearly. Therefore, the initial gravitational wave spectrum can 
correspond to those predicted by Supersymmetry, Loop Quantum Gravity, String Theory, 
Deformed Special Relativity, Variable Speed of Light theories, or many other theoretical 
models in order to provide potential tests of theoretical physics once the system is evolved 
with the turbulent matter field. 



2.2.5. Numerical Evolutions 

The development of a stable and accurate GRMHD code has been a work in progress for 
the past several years. We now have a new GRMHD code in the testing and improvement 
stage. This code is designed to perform space-time evolutions using either 2nd order Finite, 
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4th order Finite or Pseudo-Spectral differencing methods with Magnetohydrodynamic and 
pressure-less matter fields as well as various boundary and gauge conditions. The code 
can evolve the matter field independently of the spacetime in cases where a fully dynamic 
spacetime is not needed. 

We utilize the Cactus framework to develop this code. We developed an arrangement 
for Cactus that contains the GRMHD initial data, analysis and evolution thorns. This code 
handles the physics while Cactus does the 10 and Parallelization. The code is structured so 
that all the differencing is done outside of the main loops. This allows us to choose between 
several differencing techniques such as finite differencing or Fourier spectral differencing at 
run-time. We also used the Cactus Method of Lines routines to supply the time integrators. 
The spacetime (LHS of Einstein's Equations) can be ev olved using a stron gly hyperbolic 



form of the BSSN equations as defined by Brown et al (jBrown et al 



20121 ). The matter 



field (RHS of Einstein s Equations 



defined by Duez et al (IDuez et al. 



is evo lved by the form of the GRMHD equations as 



2005al ) with divergence cleaning and artificial viscosity. 
Periodic boundary conditions are also used so that the simulation domain can accurately 
represent a homogenous slice of a much larger universe. 

In addition to developing an accurate GRMHD evolution code, it is important to 
extrapolate the data so that it can be compared to cosmological observations. Gravitational 
Waves are calculated directly from the Stress- Energy Tensor using it's Quadrupole Moments. 
By doing this, the spectrum and relative amplitude of primordial gravitational waves 
created as a result of the turbulence in the matter field can be determined. Eventually, 
these results can be compared to stochastic gravitational wave data from GW observatories 
and observations of the Cosmic Microwave Background. 

There are many challenges to evolving the numerical code. First, there are the standard 
difficulties of dealing with a nonlinear code. Speed and accuracy are the most important 
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issues. Also, there are additional challenges because the GRMHD code utilizes a nonlinear 
primitive variable solver to recover elements of the Stress-Energy Tensor from the MHD 
evolution variables, shock capturing techniques and a technique called divergence cleaning 
to maintain physical values for the B-field. Optimizing and improving these solvers are 
essential to developing a fast, accurate and stable code. 



Be 



ore running the experiments we thoroughly tested the code. The Duez 



paper (IDuez et al 



2005a]) suggested four tests of a GRMHD code, however, because 
of the limited scope of this study I felt that only the following tests are necessary: 
Gravitational Wave-Induced MHD Waves, Mikowski Spacetime MHD Tests, such as shock 
tests and Consistency with the Standard Model of Cosmology. I will not include tests 
of Unmagnetized Relativistic Stars or Relativistic Bondi Flow in this paper because the 
spacetime that they are simulating lacks stars and black holes. 

This code is being developed and run on a variety of computing resources including: 
UHCL's Athena cluster, University of Houston's Xanadu cluster and University of Texas' 
Ranger cluster via the XSEDE network. 



2.2.6. Data Analysis 

The Data Analysis part of this project will focus around determining the gravitational 
wave spectrum from the simulation. A Fourier Analysis of the Quadrupole Moments of 
the Stress-Energy Tensor should yield the spectrum of gravitational waves produced by the 
turbulent matter field. Much of the data analysis work consists of the addition and fine 
tuning of new analysis routines in the code. 

Visualization of Cactus-generated data is done using a variety of open-source software 
such as Visit, xgraph, ygraph and gnuplot. Each requires implementation scripts to be 
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written. These scripts will tell the visualization program how to read the Cactus-generated 
data files. Modifications of the data include Fast Fourier Transforms (FFT) for spectral 
analysis. 

As these numerical experiments are being performed, the output is analyzed. The 
effects of variations in the density, temperature, magnetic field and initial turbulence will 
be studied in the output data. A Fourier analysis of the perturbed Quadrupole Moments 
will be performed in order to extract the spectrum of gravitational waves. This spectrum 
will then be extrapolated to give the current observable values. The result will be several 
templates of GW spectra resulting from different initial conditions. 



3. Testing the Code 



The first test that we performed involved generating Alfven and magnetosonic modes 
by gravitational waves and compar ing the results against the semi-analytic solutions from 



the Duez paper ( 



DuezetaL 



2005bl ). This semi-analytic solution is only valid for a time 



much less than the dynamical collapse time of the unperturbed fl uid. We began by using 



the same initial conditions as defined by Duez's General Example ( iDuez et al 



2005aj): 



h+(t,z) = h+o sin(kz) cos(kt), h x (t,z) = h x0 sin(kz) cos(kt), (29) 

P(0, z) = 1.29 x 1(T 9 , p (0, z) = 2.78 x 1(T 9 , (30) 

v i (0,z) = Q, S*(0, z) = (1.09,8.26, 14.4) x 1(T 5 , (31) 

h +0 = h x0 = 1.18 x 10" 4 , (32) 



where k is the wave number and we assume the fluid is unperturbed at t = 0. Results 
are shown in Figure 1. The test was performed using second order, fourth order and 
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spectral differencing as a one dimensional problem. The results shown used 200 grid points 
for the second order differencing, 50 grid points for the fourth order differencing and 32 grid 
points for the spectral differencing. Using more grid points for the fourth order and spectral 
differencing resulted in an unstable evolution. Additional tests were performed where the 
initial conditions were varied. For every variation the test proved successful, the analytic 
and numerical results proved almost identical to within a few percent. 

For our second test we utilized a FRW space-time that contained parameters based 
with those accepted for the universe at t = 10~ 6 s. The goal of the test was to determine if 
it evolved consistently with the Friedmann equations. We set the initial scale factor a = 6.6 
x 10 -14 , the initial Hubble Parameter H = 1.4 x 10 26 km/s/Mpc, the initial temperature 
T = 4.1 x 10 13 and the initial energy density po e = 1-1 x 10 23 . The system evolved until 
around t = 86s and the change in these parameters was found at the end of the simulation. 
This test again proved successful with the errors of less than one percent as shown in Table 
1. 

At this point we are prepared to add standing gravitational waves with a spec trum 



consistent to Grishchuk's predictions (IGrishchuk 



1998 



2001 



Grishchuk et al. 



20011 ). The 



gravitational waves were given random phases in order to avoid large nodes and anti-nodes. 
We then ran simulations with space-time perturbations and large (10 17 Gauss) magnetic 
fields. These simulations showed that space-time perturbations and magnetic fields had 
no significant impact on the expansion rate of the space-time and therefore adding them 
should still result in a simulation consistent with cosmological theory to within the same 
error as found in Table 1. 



Finally, we conduct ed shock tests usi ng similar initial conditions to Duez (IDuez et al. 



2005al )and Komissarov ( IKomissarov 



19971 ) as shown in Table 2. The results shown below 



are based on runs utilizing the second order finite differencing method. The fourth order 



-21 - 




Fig. 1. — Gravitational Wave-Induced MHD Waves test results. The left column used 2nd order 
finite differencing, the middle column used 4th order finite differencing and the right column used 
spectral differencing. Density is displayed here by calculating Sp/(pho), Velocity is displayed here 
by calculating 5v x I fop and B-held is displayed here by calculating 5B X / (B^ho) as in the Duez paper 



liDuez et al 



2005d ). 
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Table 1: Consistency of Standard Model of Cosmology parameters test results. 





t 


T 




H 


a 


Initial 


lir 6 


4.120 x 10 13 


1.091 x 10 23 


2.401 x 10 26 


1.0 


Final 


8.629 x 10 1 


4.446 x 10 9 


1.480 x 10 7 


2.792 x 10 18 


9.289 x 10 3 


Final/Initial 


8.629 x 10 7 


1.079 x 10~ 4 


1.357 x 10~ 16 


1.163 x 10~ 8 


9.289 x 10 3 


t a 


8.629 x 10 7 


8.586 x 10 7 


8.586 x 10 7 


8.601 x 10 7 


8.629 x 10 7 




1.077 x 10~ 4 


1.079 x 10~ 4 


1.079 x 10~ 4 


1.078 x 10~ 4 


1.077 x 10~ 4 


pe a 


1.343 x 1(T 16 


1.357 x 10~ 16 


1.357 x 10~ 16 


1.352 x 10~ 16 


1.343 x 10~ 16 




1.159 x 10~ 8 


1.165 x 10~ 8 


1.165 x 10~ 8 


1.163 x 10~ 8 


1.159 x HT 8 


a a 


9.289 x 10 3 


9.266 x 10 3 


9.266 x 10 3 


9.274 x 10 3 


9.289 x 10 3 


t b 


0.0 


4.338 x 10 5 


4.338 x 10 7 


2.824 x 10 5 


-2.980 x 10~ 7 




2.716 x 10~ 7 


0.0 


-5.833 x 10~ 16 


9.506 x 10" 8 


2.716 x 10~ 7 


pe h 


1.361 x 10~ 18 


2.933 x 10~ 27 


0.0 


4.773 x 10~ 19 


1.361 x 10~ 18 


H b 


3.804 x 10~ n 


-2.051 x 10~ n 


-2.051 x 10" 11 


0.0 


3.804 x 10~ n 


a b 


1.637 x 10- 11 


2.338 x 10 1 


2.338 x 10 1 


1.521 x 10 1 


0.0 


t c 


3.351 x lO" 3 














1.678 x 10~ 3 








pe c 






6.686 x 10~ 3 






H c 








3.357 x 10~ 3 




a c 










1.678 x 10~ 3 



a Predicted Ratio Based on the Friedman Equations using the 
column parameter 

b Difference between the Predicted and Calculated Ratio 

c Average Error: Difference / Calculated Ratio averaged among 

all nonzero columns 
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finite differencing method gives almost the exact same results when less grid points. These 
tests could not be completed with our spectral differencing method because the suggested 
shock tests are not periodic. 

Overall, the artificial viscosity methods used in this code seemed to produce less 
accurate results when comp ared with the Hi gh Resolution Sho ck Capturing Methods used 



by Duez (IDuez et al. 



2005af )and Komissarov (IKomissarov 



19971 ). Fast and Slow Shocks: The 



shock fronts for both of these cases was m ore distorted than the results of Duez (IDuez et al. 



2005al )and Komissarov ( IKomissarov 



19971 ). Also, the Fast Shock appeared to move slower 
than the Slow Shock which doesn't seem to agree with the standard results. Switch- on/ off 
Rarefaction: While the Switch-off (Fast Rarefaction) seemed to agree with the published 
results the Switch-on (Slow Rarefaction) appeared distorted at the shock front. Alfven 
Wave: Our Alfven Wave results seemed to agree with the published results except for the 
extra dip for z < 0. Shock Tubes 1 and 2: Shock Tube tests also produced distorted shock 
fronts, particularly for z > 0. Collision: The Collision test produced the best results when 



compared to the established results. 



4. Discussion 

We now have several parameters to work with in developing numerical experiments. 
Assuming the expansion properties of the background spacetime and composition of the 
matter field are fixed, we can alter the metric perturbations (scalar and tensor), initial 
magnetic field strength and the turbulence of the initial matter field. By altering these 
properties, we can perform a variety of numerical experiments to determine the effects of 
scalar perturbations and gravitational waves on structure formation, limits on primordial 
magnetic fields, properties of gravitational waves formed by a turbulent plasma, the 
dynamics of a turbulent plasma in an expanding universe and other interesting scientific 
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Table 2: Initial states for one-dimensional MHD tests. 



Test 



Left state 



Right State 



Grid tfinai 



Fast Shock u* = (25.0, 0.0, 0.0) 

(/i = 0.2 b ) SV-v/iTT = (20.0, 25.02, 0.0) 
P = 1.0, po = 1.0 



u* = (1.091, 0.3923, 0.00) n = 40 2.5 
B i /y/An= (20.0,49.0,0.0) 
P = 367.5, po = 25.48 



Slow Shock 
(p = 0.5 b ) 



u* = (1.53,0.0,0.0) 
£7>/47r= (10.0,18.28,0.0) 
P = 10.0, po = 1.0 



= (0.9571, -0.6822, 0.00) n = 200 2.0 
BVa/Itt = (10.0,14.49,0.0) 
P = 55.36, po = 3.323 



Switch-off Fast 
Rarefaction 



«* = (-2.0,0.0,0.0) 
P7V47r= (2.0,0.0,0.0) 
P = 1.0, p = 0.1 



«* = (-0.212,-0.590,0.0) n=150 1.0 
57>/47r = (2.0,4.71,0.0) 
P = 10.0, po = 0.562 



Switch-on Slow 
Rarefaction 



v* = (-0.765,-1.386,0.0) 
57747= (1.0,1.022,0.0) 
p = 0.1, po = 1.78 x 10~ 3 



v? = (0.0, 0.0, 0.0) n= 150 2.0 

57747 = (1.0,0.0,0.0) 
p = 1.0, po = 0.01 



Alfven Wave 
(p = 0.626 b ) 



u L = (0.0,0.0,0.0) 
p774¥= (3.0,3.0,0.0) 
P = 1.0, po = 1.0 



u* = (3.70, 5.76, 0.00) n = 200 2.0 
P774¥= (3.0,-6.857,0.0) 
p = 1.0, p = 1.0 



Shock Tube 1 



«* = (0.0,0.0,0.0) 

p774T= (1.0,0.0,0.0) 

P = 1000.0, po = 1.0 



^ = (0.0,0.0,0.0) n = 200 1.0 

p774T= (1.0,0.0,0.0) 

P = 1.0, p = 0.1 



Shock Tube 2 



v} = (0.0,0.0,0.0) 
P7V4^= (0.0,20.0,0.0) 
P = 30.0, po = 1.0 



u l = (0.0,0.0,0.0) n = 200 1.0 

p774^= (0.0,0.0,0.0) 

P= 1.0, p = 0.1 



Collision u l = (5.0, 0.0, 0.0) 

57747= (10.0,10.0,0.0) 
P = 1.0, p = 1.0 



u i = (-5.0, 0.0, 0.0) n = 200 1.22 

57747= (10.0,-10.0,0.0) 
P = 1.0, po = 1.0 



a In all cases, the gas satisfies the T-law EOS with T = 4/3. For 
the first 7 tests, the left state refers to x < and the right state, 
x > 0. 
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Fast Shock - Pressure 

2 -1.5 -1 -0.5 0.5 1 1.5 2 
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2 -1.5 -1 -0.5 ( 
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2 -1.5 -1 -0.5 0.5 1 1.5 2 
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""""^-rn^ Alfven Wav 


s - Pressure 


2 -1.5 -1 -0.5 ( 


) 0.5 1 1.5 2 






Shock Tube 1 - Pressure 


? -15 -1 -n 5 „„„ ( 


n 5 1 15 ? 




52-n 

^N. Shock Tube 2 - Pressure 
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2 -1.5 -1 -0.5 0.5 1 1.5 2 




Collision - Pressure 




0.5 1 1.5 2 

Fast Rarefaction - Vz 




Alfven Wave - Vz 



Shock Tube 1 - Vz 



« 1 tf. 



Shock Tube 2 - Vz 


1.01 


v v 


2 -1.5 -1 


-0.5 ( 


0.5 1 1.5 


2 




-0.99 

V 


Collision - Vz 





1.5 2 -2 



Fig. 2. — Shock Test Results shown using second order finite differencing and artihcial viscosity 
at t fi na i time. The pressure profiles are shown on the left and the z component of velocity prohles 
are shown on the right. 
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properties. 



One of the most interesting of these experiments could inv olve the interaction between 



gravitational waves and the primordial plasma. Work by Duez (IDuez et al. 



2005bJ), showed 



that gravitational waves ca n induce oscillatory modes in a plasma. According to Shebalin 



(phebafin 



2002 



2005 



20071 ). large-scale coherent structures grow naturally out of MHD 
turbulence. Here, structure is defined as strengthening magnetic fields, permanent density 
and temperature variations and secondary relic gravitational waves. One can assume that 
space-time perturbations in the early universe (sometime after t = 10 -6 secon d) interacte d 



with t he primordial plasma and resulted in Alfven and magnetosonic modes (IDuez et al. 
2005bl ). These mode s then interacted dynamica lly, possibly resulting in turbulence and 



structure formation ( iShebalin 



2002 



2005 



20071 ). Using the techniques of numerical 



relativity, we can test this assumption. 
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